The Finite Size Error in Many-body Simulations with Long-Ranged Interactions 
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We discuss the origin of the finite size error of the energy in many-body simulation of systems 
of charged particles and we propose a correction based on the random phase approximation at 
long wave lengths. The correction comes from contributions mainly determined by the organized 
collective oscillations of the interacting system. Finite size corrections, both on kinetic and potential 
energy, can be calculated within a single simulation. Results are presented for the electron gas and 
silicon. 



The accurate calculation of properties of systems con- 
taining electrons is a very active field of research. Among 
the possible numerical approaches, quantum Monte Carlo 
methods are unique in their ability to produce reliable 
ground state properties at a reasonable computational 
costQ. However, in the simulation of bulk systems, cal- 
culations are necessarily performed using a finite number 
of electrons with a consequent loss in accuracy. In or- 
der to reduce this bias, the system and, hence, the pair 
interaction are made periodic in a supercell with basis 
vectors {L a }a=i,2,3- (In the case of a crystal these vec- 
tors define a supercell of the unit cell.) This is achieved 
by using the Fourier components of the interaction at the 
reciprocal wave vectors of the supercell, i.e. G such that 
cxp(iGL Q ) = 1. Singular long-ranged potentials, such as 
the Coulomb interaction, arc computed by splitting the 
sum into a portion in real and reciprocal space 0- Al- 
though using the periodized potential reduces the finite- 
size effects, some error still remains; the one on the en- 
ergy, for example, often exceeds the statistical noise and 
other errors characteristic of quantum simulations Fi- 
nite size scaling is possible, but difficult, because the cost 
of a simulation increases rapidly with the number of par- 
ticles in the supercell. Here we present an approach that 
reduces the finite size errors. 

In Fourier space and atomic units, the electron-electron 
potential is: 
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where pg = X^i ex P(^ r i) ^ s the Fourier transform of 
the charge density, ft the volume of the supercell and N 
the total number of particles. The boundary conditions 
on the wave function can be chosen as ^(..r^ + L Q ..) = 
exp(i6 a )^(..ri..) where a is the "twist" of the phase 
in the a th direction. Periodic boundary conditions have 



6 a = 0. When there is no long range order, finite size 
errors are reduced by averaging over twists (i.e. k-point 
sampling or Brillioun zone integration) 4] . This comes at 
little cost in simulations since the average is also effective 
in reducing the statistical noise. Even when this is done, 
the expectation value of the potential energy remains ex- 
pressed as a series over G vectors and is determined by 
the static structure factor Sn(G) = (pgP-g)/N. As the 
system size increases, the mesh of G vectors gets finer 
and the series eventually converges to an integral corre- 
sponding to the exact thermodynamic limit. 

The error using a simulation box with TV particles is 
therefore given by 
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Its leading order contribution is given by the Madelung 
constant, vm-i an d corresponds to the difference 
-e 2 f(2nk)- 2 dk + 2Tre 2 n- 1 Y lG , G- 2 . It scales as l/L 
because of the omission of the G = contribution from 
the sum and its value is proportional to e 2 f v (2nk)~ 2 dk 
where I? is a domain of volume (2ir) 3 /fL 

Although vm is generally introduced using a real space 
picture, as the interaction between images, the above per- 
spective can be easily generalized to the next order cor- 
rection. The remaining part of the error is determined by 

i) the substitution of S'oo(k) by the computed Sjv(G) and 

ii) the discretization of the integral of e 2 S , (k)(47r 2 fc 2 ) -1 . 
The behavior of S(k) at large k is determined by the 
short range correlation and can be neglected. This is 
apparent if the potential is decomposed in a short and 
long range part. The long range part, whose expectation 
value is affected by the finite size, decays quickly to in 
reciprocal space so that the behavior of S(k) at large k is 
irrelevant. Moreover, in the limit k — ► 0, one knows that 
the random-phase approximation becomes exact and de- 
scribes independent density-fluctuation modes||. In the 
small k region the random-phase approximation suggests 
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S 00 (G)^S N (G) (3) 
and implies that the leading order contribution to the 
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FIG. 1: Lower panel: Static structure factor for the electron 
gas at r 3 = 10. Upper panel: AS — Sjv(fc) — 5*66 (k). The 
difference is computed using a spline function interpolation 
of See- 



error comes from point ii) above. It is an integration er- 
ror that, analogously to the Madelung constant, comes 
from the omission of the G = volume element from 
the energy sum. Scaling of the finite size errors is then 
determined to leading order by this missing contribution 
i.e. e 2 J v 5'(k)(47rfc 2 )^ 1 (ik where V is a domain centered 
on k = whose volume equals (2-7r) 3 /0. This, together 
with the characteristic quadratic behavior of S{k) for cor- 
related charged systems, leads straightforwardly to the 
well known l/fl scaling of the error. [(| Thanks to the 
validity of the random-phase approximation, 5(k) can 
be determined in the small-fc region either analytically 
or from a knowledge of the Sn (k) computed in the simu- 
lation. Once S{k) is known, one can accurately compute 
the correction. 

We looked at jcllium as a test case to judge to what 
extent Eq|3 is verified. Results for SjyQc) computed in 
variational Monte Carlo simulations at r s = 10 for 12, 24 
and 54 particles are shown in Fig^ As we increase the 
number of particles, the grid of k points for which Sn 
is defined shifts, but the values of Sn fall on a smooth 
curve, independent of N. 

Let us now consider the kinetic energy. It is impor- 
tant to distinguish between the effects due to momentum 
quantization and long range correlation. When using a 
twisted boundary condition 9 in a cubic cell, the kinetic 
energy is given in terms of the momentum distribution 

by 
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When using a single twist, for example periodic boundary 
conditions, the finite size error is, once again, composed 
of two contributions: the integration error and the error 
in approximating the exact momentum distribution, n^, 



with njy. To better understand the latter point, consider 
the fourier transform of the momentum distribution: the 
one-body density matrix. This is equal to the integral 
over particle coordinates of ^(ri + r, r2...)\P(ri, ra-..) 
and converges to the exact one as soon as the correla- 
tion length is less than the size of the simulation box. 
Under the assumptions of no long range correlation, this 
criterion is eventually met so one has n/v(k) = n^k) 
and the error comes again from approximating the ther- 
modynamic integral with a sum. At variance with the 
potential energy case, a change in the twist modifies the 
grid over which the kinetic energy is computed (see Eq^J) 
so that the error can be made arbitrarily small by in- 
creasing the density of twist angles. One can get away 
with a small number of special fc-point in the case of 
semiconductors but a finer grid is needed for a Fermi 
liquid due to the discontinuity at the Fermi surface. In 
the latter case the occupation of the single-particle states 
varies with the twist and one can use the grand-canonical 
ensemble to eliminate this source of error Q. 

Consider now the effects due to long range correlation. 
In Coulomb systems the interaction causes the wave func- 
tion to have a charge-charge correlation factor: the Jas- 
trow potential. Within the random phase approximation 
the ground state of the system is described by a collection 
of dressed particles interacting via short range forces and 
quantized coherent modes, the plasmons. Accordingly, 
the many-body wave function factorizes asQ 
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where ^ s . r . only contains short range correlations and mg 
decays quickly to as G increases and diverges as G -2 at 
small G. Because of this divergence, un converges very 
slowly to its thermodynamic value and the average over 
twists provides only a partial correction. Although one 
can address the bias on the momentum distribution^ 
directly, we here employ a different route. Thanks to 
Green's identity the kinetic energy is written as T = 
— (h 2 V 2 In ip/ 4m) 01 with a contribution coming from 
the Jastrow potential given by 
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Hence the error of the kinetic energy also has the form of 
Eq|21 a I/O finite size error in the kinetic energy corre- 
sponding to the omission of the G = term in 
This is an integration error provided does not depend 
on the system size. This must be the case whenever Eq|3 
is satisfied since a difference in Uk would necessarily im- 
ply a difference in S'(k). 

Errors in the potential and kinetic energy have there- 
fore a very similar mathematical structure. To com- 
pute the two corrections we use the Poisson sum formula 
C(L) = O -1 J2g C(G) where ( and ( are a Fourier 
transform pair. By separating the L = and G = 
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FIG. 2: Energies per particle of the electron gas at r a = 10 
in Rydberg as a function of the inverse particle number. Cir- 
cles are the Monte Carlo energies averaged over twist angles. 
Squares are the energies after the additional hui p /2N correc- 
tion (see text). 
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FIG. 3: Structure factor (left) and Jastrow potential (right) 
for diamond silicon at ambient pressure. The continuous lines 
are fit to the data (see text). The Jastrow potential shows a 
k~ 2 divergence at small k that was not explicitely imposed 
but obtained through energy variance minimization using the 
CASINO code. 



contributions from the two sums we get the expression 
for the error 

v ' J G^0 L^0 

(7) 

One sets £(U) equal to the k = limit of 2-rre 2 S(k)k 2 
or t\}k 2 u{k)/Am for the correction to the potential and 
kinetic energy respectively. 

We first apply these corrections to the electron gas for 
which the small k limits of S{k) and u{k) are known 
from the random phase approximation as, respectively, 
hk 2 /2mujp and Aire 2 / ftw p k 2 where ui p is the plasma fre- 
quency. In our tests, the wave function had a backfiow- 
Jastrow form|l2l| and simulations were performed in the 
grand-canonical ensemble. Thanks to the translational 
invariance of the Hamiltonian, the wave function factor- 
izes as exp(i6J2i r i/L)$> where the periodic part, is 
invariant in a finite pocket of /c-space around each twist 
angle. In each pocket the energy dependence on 8 is triv- 
ial and one can exploit this fact to reduce the number of 
twist angles to be the number of inequivalcnt pockets. 
This, together with cubic symmetry, drastically reduces 
the number of needed twist angles to between 20 — 200 for 
an unpolarized system with TV ~ 10 — 100. The leading 
order correction due to long range correlations to kinetic 
and potential energy are equal and sum up a total er- 
ror Ajv = hujp(2N)~ 1 . Corrected and uncorrected varia- 
tional energies are shown in Fig[3for r s = 10. Diffusion 
Monte Carlo values are uniformly shifted to lower energy 
by 0.6 mRyd/electron and show similar behavior. One 
can see that the bias due to the small size of the simu- 
lation cell is tremendously reduced, so that the TV = 12 
case is already satisfactory. 

As a second example we considered the diamond 



structure of silicon at ambient pressure (r s = 2 .0) . 
Calculations were performed using the CASINO [l3l] 
code, a Slater- Jastrow wave function, a Hartree-Fock 
pseudopotential |l4l Il5j and periodic boundary condi- 
tions. The orbitals used for the trial function (Hartree- 
Fock) were from the CRYSTAL98 code0. To eliminate 
the effects of momentum quantization we used a correc- 
tion based on the density functional eigenvalues of those 
single-particle states periodic in the simulation cell. Al- 
though this is quite common practice it involves another 
uncontrolled approximation and results depend weakly 
on the functional employed (we used the local density ap- 
proximation) . The parameters in the Jastrow potential 
and a one-body term were optimized. The two-particle 
Jastrow factor was made up by a spherical short range 
part and a plane wave expansion including 3 shells of 
k-points|l7j. One needs the plane wave expansion to 
accurately reproduce the behavior of the exact Jastrow 
factor at small k, especially in the case of small simula- 
tion cells. To further eliminate errors in the wave func- 
tion we correct the diffusion Monte Carlo value of S(k) 
by S*q MC - 5q MC which leads to an estimate correct to 
second-order in the wave function. 

For Eq0we assumed S(k) = 1— exp(— ak 2 ) and u(k) = 
4ira[k~ 2 — (k 2 + a _1 ) _1 ]|l3- When k is expressed in 
atomic units, the optimal value of a and a were found 
to be 0.72 and 1.0 respectively, leading to corrections 
of 0.13/TV and 0.092/TV hartree per electron for potential 
and kinetic energy. Results after the two corrections were 
applied are shown in Fig0] Even for the smallest cell 
(cubic, with 8 Si atoms), the error in the energy is of 
the order of 1 mHartree/electron (0.1 eV/atom) when 
compared to the value extrapolated for the infinite size. 

To conclude, we propose a way to estimate the errors 
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FIG. 4: Diffusion Monte Carlo energies per electron in di- 
amond Silicon at r a = 2.0. Energies and the S(k) and u(k) 
used to compute the correction are all obtained in simulations 
with the same number of particles. The smallest cell is the 
conventional fee cubic cell of diamond. The two intermidiate 
ones are, respectively, 2x2x2 and 3x3x3 supercells of the 
primitive cell. The largest one is a 2 x 2 x 2 supercell of the 
conventional cubic cell. 



in the potential and kinetic energy under the assump- 
tion that the low k behavior of the correlation factor 
is unchanged upon variation of the simulation cell size. 
This scheme is suggested by the random-phase approxi- 
mation that describes independent collective mode in the 
limit k — > 0. The dominant finite size errors on potential 
and kinetic energy are integration errors that can be es- 
timated by using the properties of the charge structure 
factor and the Jastrow potential at long wavelength. The 
behavior of these quantities in the small k limit can ei- 
ther be obtained analytically (e.g. for the electron gas) 
or from results with accurate optimized trial wave func- 
tions. This approach can be used to obtain energies close 
to the thermodynamic limit without performing a scal- 
ing analysis using different sized systems or assuming the 
finite-size behavior is given by Fermi liquid theory or ap- 
proximated by density functional theory. 
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by the U.S Army Research Office under DAAD19-02-1- 
0176. Computational support was provided by the Mate- 
rials Computational Center (NSF DMR-03 25939 ITR), 
the National Center for Supercomputing Applications at 
the University of Illinois at Urbana-Champaign and by 
CNRS-IDRIS. 
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